

In [1]:

    
import skyfield.api
from skyfield.sgp4lib import EarthSatellite
from skyfield.constants import AU_KM, DAY_S
from skyfield.functions import length_of
import numpy as np

Update 2020-07-17: this has been re-run with the last version of Skyfield, where the problem has been fixed (see #224). I have mantained all the comments about the problems that Skyfield had, so now they don't match the results.

As a test TLE we use a TLE for Es'hail, which is in a nearly GEO orbit at 24ºE.



In [2]:

    
sat = EarthSatellite('1 43700U 18090A   18335.89431171 +.00000133 +00000-0 +00000-0 0  9993',\
                     '2 43700 000.0858 245.4352 0001094 006.6237 164.6135 01.00274015000309')

sat.epoch.utc









    Out[2]:





CalendarTuple(year=2018, month=12, day=1, hour=21, minute=27, second=48.5317263007164)

The TLE epoch is near 2018-12-01 21 UTC. For simplicity, we compute the position and velocity of the satellite at 21:00:00 UTC and 21:00:01 UTC.



In [3]:

    
ts = skyfield.api.load.timescale()
t0 = ts.utc(2018, 12, 1, 21, 0, 0)
t1 = ts.utc(2018, 12, 1, 21, 0, 1)
rv0 = sat.at(t0)
rv1 = sat.at(t1)

The computation below should give something close to zero, as we are approximating the position at t1 by the position plus velocity at t0.



In [4]:

    
rv1.position.km - rv0.position.km - rv0.velocity.km_per_s









    Out[4]:





array([-1.08401336e-05, -1.40169930e-04,  2.14398141e-06])

The Z component should be much smaller. Indeed, let us compare the Z component of the velocity estimated by taking the position at t1 minus the position at t0, with the velocity computed at t0. We see that they differ by an order of magnitude.



In [5]:

    
rv1.position.km[2] - rv0.position.km[2]









    Out[5]:





-0.00026495564615558465



In [6]:

    
rv0.velocity.km_per_s[2]









    Out[6]:





-0.0002670996275668849

If we ask the satellite to compute its position in ITRF instead of GCRS, then the problem doesn't happen.



In [7]:

    
rv0_ITRF = sat.ITRF_position_velocity_error(t0)[:2]
r0_ITRF = rv0_ITRF[0] * AU_KM
v0_ITRF = rv0_ITRF[1] * AU_KM / DAY_S
rv1_ITRF = sat.ITRF_position_velocity_error(t1)[:2]
r1_ITRF = rv1_ITRF[0] * AU_KM
v1_ITRF = rv1_ITRF[1] * AU_KM / DAY_S

The same calculation for the Z component now gives a much smaller error.



In [8]:

    
r1_ITRF - r0_ITRF - v0_ITRF









    Out[8]:





array([ 2.56276996e-05, -6.11424304e-05,  2.16614035e-06])

If we compare the Z components of the velocity estimated by taking differences or by getting the velocity at t0, we see that they are very close.



In [9]:

    
r1_ITRF[2] - r0_ITRF[2]









    Out[9]:





-0.00453944263882633



In [10]:

    
v0_ITRF[2]









    Out[10]:





-0.00454160877917917

So it seems there is a problem with how the velocity is handled by ITRF_to_GCRS2().

Update 2019-07-24: This shows how this bug affects Doppler computations.



In [11]:

    
observer  = skyfield.api.Topos(latitude = 40.595865, longitude = -3.699069, elevation_m = 800)
los0 = (sat - observer).at(t0)
los1 = (sat - observer).at(t1)

Doppler computed as projection of velocity vector onto line of sight vector:



In [12]:

    
np.sum(los0.velocity.km_per_s * los0.position.km) / length_of(los0.position.km)









    Out[12]:





0.0005891289567179068



In [13]:

    
np.sum(los1.velocity.km_per_s * los1.position.km) / length_of(los1.position.km)









    Out[13]:





0.0005891225228701842

If instead we compute the velocity vector subtracting the positions at t1 and t0:



In [14]:

    
v = los1.position.km - los0.position.km
np.sum(v * los0.position.km) / length_of(los0.position.km)









    Out[14]:





0.00048325541071942595

There is a huge difference (a factor of 2) between both ways of computing the Doppler. The second way gives the correct result.